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Abstract. In this paper we show a simphfied optimisation approach for free boundary problems in 
arbitrary space dimensions. This approach is mainly based on an extended operator splitting which 
allows a decoupling of the domain deformation and solving the remaining partial differential equation. 
First we give a short introduction to free boundary problems and the problems occurring in optimisation. 
Then we introduce the extended operator splitting and apply it to a general minimisation subject to 
a time-dependent scalar-valued partial differential equation. This yields a time-discretised optimisation 
problem which allows us a quite simple application of adjoint-based optimisation methods. Finally, we 
- - - verify this approach numerically by the optimisation of a flow problem (Navier-Stokes equation) and the 

final shape of a Stefan-type problem. 
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1. Introduction 



Optimisation of free surface problems [5J I16j often occurs in industrial applications. Some examples 
are the stabilisation of a liquid surface for sloshing [TT] or optimising the shape of solidification processes 
[10) . Free surface problems are still challenging from an analytical as well as a numerical point of view. 
Here, the domain is an unknown of the equation system which depends on the states, e.g. a water surface 
is driven by the flow velocity Since these problem arc already hard to handle, the optimisation of such 



(~| ' processes is very complex. Especially adjoint-based approaches [18] arc very difficult to apply due to 



the statc-dcpcndcnt domain. Here, several assumptions and methods were derived to handle this kind of 
problem. For special cases it is possible to describe the free boundary by a graph [15] or introducing a 
level-set or phase field function [2]. Another approach is the puUback of the time- and statc-dcpcndcnt 
domain to a reference domain and perform all calculations in there. From an optimisation point of view, 
all of these methods have the disadvantage of very complex derivatives describing the variation of the 
^ ' domain. Note that often these derivatives are, in contrast to stationary problems, hard to interpret for 

time-dependent problems. 

' In this paper we show a simplified optimisation approach for free boundary problems in arbitrary space 

CO . dimensions which bases on an extended time-discretisation of the problem. We consider the problem: 

' Minimise J{y,il,u) subject to the free surface problem of finding (y, O) such that 

, dty{t) + A{t)y{t) = f{u{t)) in n{t) y{0) = yo in 0(0) 

^' B{t)y{t) ^ g{u{t)) onr(i) ^(y, O) = 

holds for all t G (0,r]. Here, 0(i) C denotes the time- and state-dependent domain, y a scalar- valued 
• function, A an arbitrary differential operator and C(y, O) a constraint function defining the free boundary. 

I Moreover, B and p : R — > R denote appropriate boundary conditions, u a control function and / : K — s- M 

■ a right hand side term depending on u, e.g. a localisation function. The difficulty is the dependency of the 

domain 0(t) on the solution y of the partial differential equation. To obtain the domain, we solve, roughly 
speaking, a minimisation problem for each time step t G (0,T] in order to fulfil the constraint for the free 
boundary. Hence, the minimisation of the cost functional J would be subjected to the minimisation of 
the constraint function C in order to solve the state equation ([T|). 

In order to avoid the minimisation problem for solving the state equation, we perform a complete 
puUback of equation ([1]) to a reference domain O. For this, we introduce a flux function : K ^ K"^ 
depending on the state y. This function is given by a characteristic velocity (e.g. flow) or a smooth 
continuation of the boundary motion, cf . figure [1] Together with the transformation $ given by 

dt^ = F{yo^) inOx(0,T) with $(A,0) = X in O 

we resolve the constraint function for the domain by 0(i) = $(0,t) and hence C{y, $(0)) = holds, see 
figure [21 Therefore, we reformulate the original minimisation problem of the cost function J to: Minimise 
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fixed boundary 



Figure 1. Artificial convec- 
tion. The arrows illustrate the 
flux function F which provides 
the deformation field to gener- 
ate the new domain or transfor- 
mation $. 




Figure 2. Fullback to a refer- 
ence domain by the transforma- 
tion $. Capital letters, e.g. X 
denote coordinates in the ref- 
erence domain (left) and lower 
case letters, e.g. x, denote coor- 
dinates in the transformed do- 
main. 



j{y, u) subject to 

dty + A{^)y = f{uo^) in £1 x (0,T) 
(2) B{^)y = g{uo^) in f x (0, T) 

y{0) = Vq in 

with y :=y o $ and lA^{x) x as the identity map in (l. Note that also the operators A and B change 
to and -B($), respectively. Now the domain is fixed but the dependency of the differential operator 

on $ yields very complex derivatives. 

In the following we derive a simplified optimisation approach for free boundary problems using an 
extended time discretisation which is based on an operator splitting. This approach is a blend of the 
original ([T|) and transformed formulation ([2]) of the optimisation problem. 

2. Formulation of the State Equation 

In this section we reformulate the state equation by a time-discretisation. Particularly this is done by 
applying an operator splitting scheme, which allows a simplified optimisation approach later on. First, 
we show a first order splitting scheme for a simple convection-diffusion equation in and transfer the 
results to the state equation ([T]). Finally, we introduce an appropriate time-discrete Hilbert space which 
simplifies the application of adjoint-based optimisation. 

2.1. Basics of Operator Splitting. In the following we use a Yanenko splitting [7] which is illustrated 
by a simple convection-diffusion equation in 

dty + V ■ Vy - Ay = f in R'* x (0, T] with y{0) = yo in R'^ 

where y : M'' — !■ M and i; : M'' — > denotes a smooth vector field. Note that the velocity field can also 
depend on y. This problem is divided into two subproblems. For the first time interval [0,t] the Yanenko 
splitting reads 

dty* + v • Vy* = in R'^ x (0, r] dty - Ay ^ f inR'^ x (0, r] 

y*(0) = yo inK'' ^ y{0) = y*{T) inR^ 

The resulting subproblems can be solved by different methods. In particular, wc consider the convection 
part from a Lagrangian viewpoint, cf. [S], which yields 

dty* =0 in x (0, r] 9*$ = vo^ in x (0, r] 

y*(0) = yo inM'^ $(X,0)=X in M'^ 

for y* :— y*o$. Thus the solution y* at time r is given by 

y*mX,T),T)^y*iX,T)^yo{X) 



dt^ = F{y) infix(0,T) 
<I'(0)=Id^ inn 
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and hence ?/*(t) — yo{^ ^ {''')) where $ ^ denotes the inverse of The diffusion part is written as 
dty-Ay = f in $(M^t) X (0,t] 

2/(0) = yo(3'"'(T)) in $(M^t) 

Note that $(M'',r) = M** but the metric, needed for example for spatial operators, is induced by the 
transformation Applying this scheme iteratively, we solve the convection-diffusion equation for all time 
intcrvalls. Finally, a time-discretisation is performed. Since the Yanenko splitting is of order 0{t), a first 
order scheme is sufficient. 

i(y("+i) _ j;(")o($("))-i) - A2;("+i) = in 

with 0^"+^) = ^l^") This approach, i.e. solving the convection part by a Lagrangian viewpoint, is 
the main principle of particle methods, cf. [13) . which will be used for the numerical results later on. 

Remark 1. Due to the definition, the continuous transformation $, given by is a diffeomorphism, cf. 
[14| . For small time steps also the time-discrete transformations (f>'-"^ are diffeomorphisms. This implies 
that det(V<f>'-"') > is always satisfied. 

2.2. Time Discretisation of the State Equation. Now wc apply the above splitting scheme to the 
state equation ((T|). For this, we extend (HJ by a convection part given by a smooth continuation of the 
boundary motion. Particularly, we obtain 

dty{t) + Fiyit)) ■ Vy{t) - F{y{t)) ■ Vy{t) + = f{u{t)) in n{t) 

B{t)y{t) - g{u{t)) on dn{t) 

y(0) = yo in n{0) 

for all t E (0,T]. Here the flux function F : M — )• M'' defines the deformation of the entire domain, cf. 
equation ([2]) and figure [T] Applying the splitting scheme of the previous section to this equation we obtain 
the time-discrete system 

(3) 



_y(")o($("))-i) - Vy^"))o($("))-i-|-.4(t("+i))y("+i) =/(«("+!)) inl]("+i) 

= 5(u("+i)) on 

with $7^"+-'^) := r2(t("+^)) = Note that the transformation, which generates the new domain, 

is solved explicitly. For this reason it is sufficient to treat the artificial convection term F{y) ■ \7y also 
explicitly in the above equation. Moreover, the transformation depends on y^"^ only. Hence we treat 
as function of y*^"^ given by 

(4) $("){2/}:-Id,,(„, +rF(y(")) 

in the following. To consider weak formulations later on, wc introduce the function space 

(5) ^ - n ^^"^ 

where F^"-* denotes the spatial space at time i*-"-*, for instance V^*^"^ := H^{fl^'^'>). These spaces are built 
recursively, that is, 

^(0) ^ -^(1) ^ -^/(2) 
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and depend on the data, e.g. right hand side or boundary conditions. If y*-"^ are Hilbert spaces, we define 
the corresponding inner product by 

(6) (x,j/)v:=E^(^^"''2/^"^)v(") 

for x,y £ V. Note that for small time steps every transformation is a diffeomorphism which maps from 
as state in the above remark. Therefore, the domain does not become singular and hence 
if is a Hilbert space then also is a Hilbert space. Consequently, © is an inner product and 

hence also F is a Hilbert space. 

3. Optimal Control Problem 

Let the state space V and the space of Lagrange multipliers Z be Hilbert spaces which are defined 
analogous to ([5]). Let the space of control functions J7 be a Hilbert space. Moreover, we introduce the 
space 

Nt 

H Yl ^^"■^ ^it^ •= L'^i^^"^)- 

n=0 

The index F denotes the corresponding function space on the boundary, e.g. iJp"'' = L^(9n'^")). 

3.1. Weak Formulation. We generalise equation (|3]) in a weak sense using the function spaces defined 
above as 

(7) ^(y^"+^^ - 2/^"^o($("){y})-l) - (i^(y("^) • V2/("^)o($("){y})-l + = b'-''+^^u 

in The operator : F^") (\/("))* denotes the weak counterpart of A and B used in 

equation ^ and B'"-' G jO.{U; of / and g. For example, we obtain for the Laplacian A^")?/*^"' :— 

Ay^"-* with Neumann boundary, f{u) ~ and g{u) = u 

(A("'y'"\A)(^,(„,).,^(„, :=- y" V2/(")-VAda; and (^("'li, A)(vW).,yw := J uXduj{x) 

nM af2(") 

for A e V^("). 

3.2. Minimisation Problem. We consider the minimisation problem 
(■p) min J{y,u) subject to e{y,u)=0 

(y,u)eVxU 

where J : V x U is a cost functional and e : V x U — !> Z* is determined by ([7]), i.e. 

Nt 

{eiy,u),X)z'.z := ^ [(2/^"^ A("))^(„, + T(A(")y("), A("))(y(„,).,^(„, - t(S(")w, A("))(^(„,).,^(„, 

n=l 

- E [(/"^°(*^"H2/})"',A^"+^))^(.+i) +r((i^(2/(")).Vy("))o($("){2;})-i,A("+i))^(„+i, 

Tl = 

Note that the constraint function e does not depend on the domains il^") explicitly. The condition 
^(n+i) _ needed to establish the space is given implicitly by the definition of 

Remark 2. On the one hand, the missing condition ~ (fi'"') in the above contraint function 

yields an underdetermined system. On the other hand, the solution of the approach introduced in section 
\2.2\ satisfies the constraint function. This solution is used to establish the optimality condition and hence 
we assume the domains to be known. In standard approaches, fixing the domain would yield no information 
about the variation of the domain or the quantity discribing it, e.g. a tranformation. Here, this information 
is not fully lost as we still obtain information about domain variations by the push forward terms in the 
constraint and cost function. 
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The above minimisation problem is the time-discrete counterpart to the minimisation of J(y, f2, u) 
stated in the very beginning. Here, aU spatial dependencies, e.g. integration domains or evaluation posi- 
tions, in the cost functional are replaced by the corresponding transformation (f>'^"'{y}. More details about 
the reformulation of the cost functional can be found in section H] or [T3] . 

To find a minimum of (V) we use the Lagrangian multiplier theorem, that is, we determine the critical 
points of the Lagrange functional L: VxUxZ-^M. defined by 

L{y,u,X) := J{y,u) + {e{y,u), X) z* .z 

Then the Karush-Kuhn- Tucker system reads 

{dyL{y,u,X),-ipy)v,v ^ 0, {duL{y,u,X),ipu)u',u = and {dxL[y,u,X),^x)z*,z^O 

for all i/jy G y, G [/ and iIjx ^ Z. For more detail we refer to [18]. The partial derivatives in V and Z 
are interpreted as 

Nt 
ri=0 

The derivatives with respect to the control u and the states y*-"-* are straight forward except for the push 
forward terms, whose derivatives are derived in the following. 

Lemma 3.1. Let C M'' and <I> : 51 — S- M'' be a smooth diffeomorphism. Then for d ~ 1, 2, 3 

div(det(V$)V$"^) = 

holds. 

Proof. See [B], p. 27. □ 

Lemma 3.2. Let Vl dW^ and ^ be a smooth diffeomorphism. Moreover, let y : ^ R and X : $(0) M 
be sufficiently smooth. Then 

y(Ao$)det(V$)rfa;)[V'] = J y (Ao$) det(V$)V$"^n • V' da;(a;) 

- j det(V$)V$"^Vy Ao$ • -0 da; 

o 

holds for i^:VL^W^. 
Proof. 

d^(^JyXo<Pdet{V<^)dx^ [ip] ^ J y (£)Ao$)V;det(V$) + y Ao$det(V$)V$"^ : V^/'da; 
n n 

J y (£)Ao$)V'det(V$) - det(V<E>)Vi>"^Vy Ao$ • - ydiv (det(Vi>)V$"^) Ao$ • 

det(V$)V$"'^V$'^(VAo$)-0da; + J y(Ao$) dct(V$)V$"'^n • V dw(a;) 

dn 

by using integration by parts. Applying lemma [3Tl we obtain 

= J y(Ao$)det(V$)V$"^n--!/'da;(a;) -y det(V$)V$"^Vy Ao$ ■ da; 

on n 

by using (DXo^)^' = (VAo<l>) • i/' and div(ai3) = adiv(i?) + B{Va) for a scalar-valued function a and a 
matrix- valued function B. 

□ 



n 

- y 
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Theorem 3.3. Let ip : fi^") — > M and X : r2("+i) — > M &e sufficiently smooth. Moreover, F : y ^ F{y) 
denotes a vector field depending on y*-"', ^'-"'{j/} is given by @ and Jl'""''^) ~ Then 



D 

y ■ ■ 



'^(„,( j ((po($("){2;})-i)Adrr)[V] = j r(^(Ao$("){y})Z?^^(2/("))7A.nda; 

(Ao$("){?/})V^ ■ DF{y^'''>)i:duj{x) + ©(t^) 



/loZds /or V : ^ M. 

Proof. Due to definition we define 

((po($("){y})-i)Ada; = y (p(Ao$("){j/})dct(V$("){y})da; A:($("){y}) 
The chain rule yields the variation with respect to y'") as 

Using the definition of yields 

and applying lemma 13.21 gives 

5^(„,/C[V'g = J ^(Ao$("){y}) det(V$^"){j/})V$("){y}-^n ■ TDF{y^"^)^y dx 

det(V$)V$~^V(p(Ao$("){2/}) ■ TDF(y^'^'>)iPyduj(x) 

Since the determinant and inverse of a matrix satisfies 

det(/ + £A) = l+eTr(A) + 0(e^) and {I + eA)-^ = I - eA + Oie"^), 
respectively, for small £ > 0, we obtain 

det(V*("Hy}) = det(/ + rF(y("))) = 1 + r div + Oir^) 

and 

(V$("'{y})-i = (/ + rVi^(2/(")))-i = / - tVF(j/(")) + Oir^) 
for small time steps r > 0. Therefore, we get 

a,(„,/C[^,] ^ J ^(Ao$('%})(l + rdiv(F(y("))))(/-TVF(j/("'))'^n.r^F(y("')V',cia; 

(1 + r div(i^(2/(")))) (/ - rVF(y("'))^V(p(Ao$("){2/}) • TDF{y^''^)ijy doj{x) + 0{t^) 
which finally yields 

5^(„,/C[^,,] = y r^{\o^^^\y})DF{y^'^^)^yndx 

-(Ao$("){j/})V(p • DF{y('''>)i'y duj{x) + ©(r^) 



□ 



Remark 3. The above theorem holds for arbitrary flux functions F, i.e. F : — > W for an appropriate 
Banach space W . 
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The following corollary shows the application of the general theorems stated above to problems given 

by 

Corollary 3.4. Let y e V , ip £ F^") and X G V^"\ Moreover, let F : R ^ R'' and <^^"\y} be given by 
(HI). Then 

(9^(V3o($("){?;})~l)[?/v], A)^(„ + i) = (Ao<i>("){y}>^)^,„, 

+ r(Aoci>(»){y}F'(2;(")) • V2/("),7A^)^(„, +0(t2) 

holds for Tpip g y*^"' and 

(9,(„,(y.o(cf>(«){y})-i)[^g,A)^(.+i, =r(^(Ao$("){y})F'(y(")).n,V',)^(., 

- T((Ao<i>("){y})V^ • ^^'(2/(")), + 0{t^) 

for iPy G 

Proof. The first part is given by 



V'^(Ao$("){y}) det(V$("Hy}) = y + Tdiv(F(y(")))) da; + ©(r^) 

n<") a(-) 
which yields the assumption with div(i^(?/("))) = F'{y'^"-^) ■ \7y^"\ 



The second part is a direct consequence of theorem 13.31 bv using the fact that F : M — > M'*. 

□ 

3.3. Adjoint System. Wc apply the above theorems to the optimal control problem (V) . The derivative 
with respect to the control function u is given by 

{duL,i^u)u',u = {duJiy.u),4>u)u'M - ^r((B("))*A("),^„)c/.,c/ 

n=l 

for all ipu G U. The first variation of the Lagrange functional with respect to y^"-* is 
(9j^(„)i,V'y)(yM).,y(") = J(y,w)>j,)(y(„)). y(„) + (A("\Vy)//(") +r(i:>A(y("))V'y,A("))(y(„)).^y(„) 

+ r(F'(y(") • Vy(")(A("+i)o$("){y}),^,)^<„,] - {r(y(") (A("+i) o$H{y})F'(y(")) • n, ^,)^(„, 

- r(F'(y(")) • Vy(")(A("+i)o$("){y}), ^,,)^,„, } + ©(r^) 

for all tpy G T/^'") in the corresponding time step n = 1 . . . A^t — 1 by using corollarv 13.41 Only the inner 
products in the {} brackets arc a result of the implicit variation with respect to the domain. Terms, which 
are handled explicitly in the time-discretisation, e.g. TF(y) ■ Vy, have variations of order which are 
directly included in O(r^). All remaining terms are due to the variation of the partial differential equation 
as usual. We simplify the above result for dy{„)L as 

+ r(i?A(y("))V„A("))(^(„,).,,.(„, - (A("+l)o<i>(«){y},V,)j,(„, 

- r(div(V',F(y("))), A(" + l'o$("){y})^<„, - r(y(") (A'^+i' o$("){y})F'(2/(")) • n,^,)^(„, + ^(r^) 
Furthermore, wc obtain for the final time step Nt 

{dy{Nt)L,ll}y)^y{Nt)y y{Nt) = {dy{Nt)J{y, U),ll;y) f^y(Nt)y yiNt) _H-(Nt) 

+ r(i^A(y(^')M„A(^')) (V(Nt))._y(Nt) 
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Forward problem ► 




■• Backward problem 



Figure 3. Basic principle of the time-discrete optimisation of free boundary problems. The 
state equation yields the sequence of domains f2'"' as a consequence of $^"\y}. The adjoint 
system goes backwards, starting from fi*^'^ to £7'*^'. 

All terms of order r are a consequence of the time-implicit scheme we chose for the discretisation and 
they disappear in explicit schemes. Since we only consider small time steps r <C 1, we neglect them in the 
following for the last time step. Using the fact that 

(div(V',i^(2/("))),A)^(„, = (A^^(y("))-n,V,)^(., - (VA • F(y("))>,)^(„, 

hold, we can, roughly, identify the adjoint equation as 

_ A("+i)o$("){y}) + Z?A(y("))*A(") + -a„(„, J(y,w) = V(A("+i) o$("){y}) . F(y(")) 

T T " 

in (\/("))* and hence an equation in fl'-'^l Here C^") : ^ is given by 

(C(«)(A("+i)ocf>(«){j/}),V-)(v(.)).,v(., :=((A("+i)ocf>(«){y})(y(«)F'(2/("))+^^(2/("))).n,V')ff(.) 

and represents the boundary values. It is, among others, a consequence of the implicit domain variation. 
The gradient of the reduced cost functional J{u) :— J{y{u), u) is given by 

VJ(«) := 7^^l (^^duJ{y,u) - X;(B("))*A(")^ 
where TZjj denotes the Riesz isomorphism TZjj : U ^ U*. 

3.4. Conclusion. The optimisation approach described in this section is very close to the numerical 
implementation of free boundary problems. The basic principle is illustrated in figure [31 Here, the state 
(or forward) equation yields the domains fi'-"-' recursively by the transformation $'^"^{y}. Consequently, 
the adjoint (or backward) equation uses the same domains, starting from fi'^') to n'^^\ Particularly, the 
push forward term of the discrete time derivative in the forward problem changes to a puUback term in 
the adjoint equation. 

Instead of deriving the variation of each domain il'") directly, it is handled by the variation of the 
push forward terms, e.g. in the time-derivative, with respect to $^"^{y} and hence to y. This procedure 
allows a easy derivation of the adjoint equation as only a few simple (explicit) terms are needed to obtain 
information about the shape dependency. Moreover, we showed above that terms, handled by an explicit 
time-discretisation, do not yield a contribution to the domain variation as the derivatives are mainly of 
order t^. The next section shows some numerical examples applying the above method. 

4. Numerical Examples 

In this section we verify the optimisation approach shown in the previous sections numerically. The 
first test case involve the Navier-Stokes equation with a free surface, c.f. [T], where the transformation is 
given in a natural way, that is, the flow velocity is used. The second example is a Stefan-type problem 
where the final shape of a melting or solidification process is optimised. Here, the transformation is not 
given by a flow velocity but an artificial motion, similar to an ALE [3] method. 
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4.1. Navier-Stokes equation. In this example, we optimise the filhng of an open hquid tank with small 
obstacles shown in figure 21 In particular, we search an inflow profile such that for a given outflow profile 
the free surface remains as still as possible. This yields the optimisation problem: Minimise 

T 

J{u,n) :=^y j \u[x,t)\^ du{x)dt 
r* 

subject to the Navier-Stokes equation 

(9a) dtU + u-Vu- i^Au = - Vp in fit x (0, T) 

(9b) divM = infitx(0,r) 

(9c) u = onT^x{0,T) 

(9d) u = Uo on To x (0, T) 

(9e) Vp-n^c onrjX(0,T) 

(9f) VM-n = onr}x(0,r) 

(9g) p = onr}x(0,r) 

(9h) m(0) = in no. 

The domain Vlt depends on the velocity profile u at the free surface F^. In particular, the motion of 
is given by u\-pt . For convenience we use (j9cl) as inflow condition. Here, c denotes the control variable 
depending on space and time, i.e. c : x [0, T] M. From a physical point of view, Vp ■ n is proportional 
to the inflow flux in normal direction. The condition for the free surface, (|9f[) and ( |9gP , is described in 
more detail in e.g. [8]. 

4.1.1. Time Discretisation. For the time-discretisation a Chorin projection [H [H] is used. In particular, 
we choose the domain transformation ^'^"^ : f]*^"^ — > as 

:=Ido(") +rM("). 

Then, the time discretisation of © reads 







in 




= idiv 


in 




= 


on F^ 




= Uo 


on To 


pin+l) 


= 


on f("+i) 


Vp("+^) ■ n 




on Fj 







in 
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where ct > denotes a regularisation parameter. All boundaries not mentioned above are set to uniform 
Neumann condition, i.e. Vp • n = or Vit ■ n ~ 0, which is a consequence of Chorin's projection. The 
cost function is discretised straight forward 

J{u,c):=-Y, r / \u^-^\^doj{x). 

4.1.2. Adjoint system. Using the previous optimisation approach, we obtain the identification of the ad- 
joint system as 



-(aL")-(aL"+^^o$("){w})) 





= VA*,") - (A(:'+i)o$("){m}) . 




in 




^divAL") 




in 




= - (l6(") • (A(,"+1'o$(" 


\u]))n 


on i ^ 




= 




on rj."^ 


z/VA^") n 


= -A(,")n -(A^I'+i^o^ 


("){w}))n 


on r, 


cVAp • n 


= -AL"' • n 




on r, 


Xin) 


= 




on Tu,,o 


aVXp ■ n 


= 




on Tu,,o 


XiN.) 


= 




in 0(^') 



The gradient of the reduced cost functional J(c) := J(m(c),c) is identified by 

VJ(c) -aApIr,. 

A detailed derivation of the adjoint system can be found in [T3]. Note that only the puUback terms on 
the right hand side, e.g. (A^I'+^^ o $("){«}) ■ Vtt^"), are a result of the domain variation. 

4.1.3. Numerical Results. For the numerical implementation we use a meshless particle method, cf. |13| . 
These methods arc superior to mesh-based methods like finite elements as no connectivity information 
is used and hence no expensive remeshing is needed if the positions of the supporting points change. 
Particularly, we use the finite pointset method, cf. |17j . The optimisation process is performed by a 
second order BFGS method with Armijo rule, cf. [S]. 

The final time is set to T = 2.5, the step size r = 0.005. The domain has a width and height of 5.0. 
Moreover, the viscosity is set to = 10 and a = 0.005. The outflow velocity is given by a parabolic profile 
with Umax = 3.0 in order to avoid singularities at the corners. 

Figure [6] and [7] show the results for the uncontrolled and optimised case, respectively. Here the colour 
represents the velocity magnitude. The uncontrolled case, that is, no inflow is given, yields a high velocity 
at the free boundary and therefore a large deformation of it, as expected. The resulting behaviour of 
the free surface is an effect of the high viscosity of the fluid. The controlled case forms out a straight 
flow, as good as possible for the given geometry, from the inflow to the outflow, which docs not affect 
the free surface strongly. Hence, the surface velocity is much smaller than for the uncontrolled case, cf. 
figure m Note that due to the setting, the desired value of the velocity it = at the free surface is not 
reachable and hence min J = cannot be expected. The corresponding evaluation of the cost functional 
and gradient norm is shown in figure [5] The gradient norm shows a strong decrease in the first iterations, 
then it becomes flatter with order 1.04 to 1.25. Note that this is not surprisingly as the problem is highly 
non convex. On the one hand, the Navier-Stokes equation yields a complicated constraint, on the other 
hand, the adaptation can cause problems for the optimisation process. 
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4.2. Stefan-type Problem. The second example deals with a Stefan-type problem ^9\. These problems 
describe phase translations, e.g from liquid to solid, see figure [5] Particularly, we solve the heat equation 

dti9 - Ai9 = / inntx (0, T) 

^ = ^m onrtx(0,r) 

^ ' i?(0) = i9o in ilo 

v„ = py-d ■ n 

with e M as the boundary temperature (melting point) and i!)o as the initial temperature distribution. 
Vn denotes the velocity of the boundary in normal direction and hence also defines ilf Moreover, /3 G M 
denotes a material constant. Note that the sign of /? indicates whether ilt is a liquid or solid phase. 

For the minimisation we start with an arbitrary domain VIq with boundary Fg and want to achieve 
a desired shape with boundary F^ at final time T. In the following, we choose an ellipse for F^. This 
desired boundary shape is parametrised in order to formulate an optimisation problem. The easiest way 
of a parametrisation is the implicit definition of an ellipse given by 

in two dimensions. Here, a and b denote the stretching in the x- and y-direction, respectively. Hence, the 
domain Qt has the boundary F^ if (jlip is satisfied for all points on the boundary F^. The minimisation 
problem can therefore be interpreted as: Find a heat source / such that 

(12) 

Tt 

is minimal subject to (fTU|) . The integral is minimal, in particular zero, if all boundary points satisfy ([TT|) . 
For all other settings, i.e. a measurable amount of points do not lay on the desired boundary, the integral 
is greater than zero. 



4.2.1. Time-discretisation. For convenience, we replace the Dirichlet boundary condition of equation (jlO[) 
by a Robin condition, i.e. 

■d = -dm 1? + kVz? ■ n = 

on Ft where k > denotes a small regularisation parameter depending on the discretisation. Since we 
are dealing with a problem without convection and we only know the motion of the boundary, the choice 
of the transformation is not obvious. In particular, we choose a smooth continuation of the boundary 
velocity by solving a Laplace equation similar to the ALE method, cf. [3|. 

—Am = in fit 

u + kVu ■ n = P{Vi} ■ n)n on F/ 

where we use a regularisation of the Dirichlet boundary condition as before. 
We obtain the time discrete system by first solving the velocity equation 

-Am(") = in f|(") 

4- kVm(") • n = ^ - )n on F^") 

then the transformation 

=Idfj(„, -|-tm(") infi(") 

and finally the heat equation 

_ (i5(")o($(«))-i)) _ A?9("+i' = • Vi?("')o($("))-i + xc("+^^ in 

T 

+ kVi?("+^) ■ n = on F^^+i' 

with the initial condition 
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where x denotes a spatial localisation function depending on x only. Note that we use the Robin condition 
of the heat equation to replace Vi? ■ n in the velocity equation and add a convection term to the right 
hand side of the heat equation due to the artificial transformation of the domain. 
The discrete cost functional reads 

J(-it,c):=i J (($(^'){it})^E($(^'){M}) 



with 



E 



1-2 

6-2 



which corresponds to Note that we use the power of 2a in the cost functional to be able to consider 
also LP norms. 

4.2.2. Adjoint system. Using the previous optimisation approach, we identify the adjoint equation and 
gradient as 

A(7' + kVA^T^ ■n = on T^") 

^ (A^"^ - A^"+^^ o$("){m}) - AA^T^ = -u^"^ ■ V(A^"+^) o$("){m}) in 

A*;') + «VA^"' • n = Ac(A^"+'' o$("){m})m(") • n - ^A(,") • n on T^") 

ioi n = 1, . . . , Nt — 1. For n = Nt we get 

-AA(f = in 

A(f') + KVA(f*) • n = -2aK(($(^'){tt})^E (^^^'H-"}) - l) E on T*^') 

and 

A*,^') - tAA*,^*) = in 



^^^^ Ar'^ + -VAf*).n=-^A(f*) 

The gradient of the reduced cost functional reads 



on 



Again, a more detailed derivation of the adjoint equations can be found in |13j . 

4.2.3. Numerical Results. The spatial discretisation is, similar to the previous test case, performed by a 
meshlcss method with adaptation. Again, the BFGS method with Armijo rule is used for the optimisation. 
The setting is a = 2, /3 = —1, k, — 0.01 and dm = 0. Moreover, we set the time step size r — 0.01 and the 
finial time to T = 0.3. The localisation function x is given by 

Ny 

x{x) c(") := E E i"^ ( - 25((^ - ^^jf + iy- y^,f] 

i=l j=l 

for the supporting points —1, —0.6, —0.2, 0.2, 0.6, 1.0, i.e. Xn = —l,yii = —1 till xee = l,j/66 = 1- 
Furthermore, the initial domain is given by a square with edge length 0.6 and the desired shape is given 
by a circle of radius 1 which yields a = 1.0 and b = 1.0. Note that this example does not base on a 
physical setting as we only want to show the feasibility of our method. 

The results are shown in figure [TUl where the colour represents the absolute value of the source induced 
by the control function, i.e. xc. Moreover, the black line denotes the desired shape at final time. The 
optimised case shows a small deformation in the first time steps, which becomes larger close to the final 
time. The shape at final time is very close to the desired one. The convergence, shown in figure [TTl is 
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approximately of order 1.5 and no Armijo step size reduction is needed. Furthermore, the plateaus in the 
the cost functional and gradient norm are due to a reset of the BFGS matrix, implemented for stability 
reasons. Note that we neglect these plateaus for the consideration of the convergence order. 

5. Conclusion 

In this paper we presented a simple approach for the optimisation of free boundary problems which is 
very close the their numerical implementation. In particular, we applied an extended time-discretisation 
based on an operator splitting. For this a decoupling of the deformation of the domain and the solution of 
the remaining partial differential equation was performed. With this approach, adjoint-based optimisation 
can be applied easily. Since we do not perform an explicit derivative with respect to the domain, this 
information is partly obtained by the variation of the push forward terms needed for the discrete time 
derivative. The numerical results based on the Navier-Stokes equation and a Stefan-type problem showed 
that this optimisation approach yields good results and is very easy to implement numerically. 
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Figure 8. Cost functional and 
relative gradient norm. 



Figure 9. Mean surface velocity 
||ix||£2(p^.) over time for the uncon- 
trolled and controlled case. 
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Figure 10. Optimal solution for the Stefan problem with a = 1.0, b = 1.0. The time steps are 
equally spaces from t = to T = 0.3. The colour represents the magnitude of the source and 
the black line the desired shape at final time. 
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(a) Cost functional. (b) Relative gradient norm. 



Figure 11. Cost functional and gradient norm for a = 1.0, b = 1.0. 



